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ABSTRACT 

The new 2-planetary system around HD 169830 has been announced during the XlX-th IAP 
Colloquium "Extrasolar Planets: Today & Tomorrow" (Paris, June 30 - July 4, 2003) by the 
Geneva Extrasolar Planet Search team. We study the orbital dynamics of this system in the 
framework of the -/V-body problem. The analysis of its orbital stability is performed using the 
long-term integrations and the fast indicators, the Mean Exponential Growth factor of Nearby 
Orbits and the Frequency Map Analysis. The HD 169830 appears to be located in a wide stable 
region of the phase space. The ratio of the mean motions of the planets HD 169830 b and c 
is between low-order mean motion resonances, 9:1 and 10:1. The long-term integration of the 
coplanar configurations, conducted over 1 Gyr, reveals that the eccentricities of the companions 
vary with a large amplitude ~ 0.4 — 0.5 but there is no sign of instability. The orbital parameters 
of the planets resemble those of another 2-planetary system, around HD 12661. Both of them 
can be classified as hierarchical planetary systems. We investigate whether these two exosystems 
are dynamically similar. Such similarities may be important for finding out if the formation and 
subsequent orbital evolution of cxoplanetary systems obey common rules. 

Subject headings: celestial mechanics, stellar dynamics — methods: numerical, N-body simulations — 
planetary systems — stars: individual (HD 169830, HD 12661) 

1. Introduction 

The radial velocity survey for extrasolar planetary systems conducted by the Geneva Extrasolar Planet 
Search Team has recently revealed a new planetary system around HD 169830. This discovery has been 
announced during the XlX-th IAP Colloquium "Extrasolar Planets: Today & Tomorrow", (Paris, June 30 - 
July 4, 2003). In this paper we perform a preliminary analysis of the orbital stability of the HD 169830 system 
and its dependence on the orbital parameters that are unconstrained by radial velocity (RV) observations. 
In the numerical experiments, we use two sets of the orbital parameters from the 2-Keplerian fits to the RV 
measurements. The first set was published by the authors on their web site in June 2003 3 and the second, 
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most recent set is from Mayor et al. (2003). Our early analysis of the first set of the orbital elements revealed 
a curious, close similarity to the orbital elements of the HD 12661 system (Fischer et al. 2003; Gozdziewski & 
Maciejewski 2003). These two hierarchical planetary systems (HPS) appeared to have very similar ratios of 
the minimal masses \i — m^/mc ~ 1.3, the semi-major axes a = ab/a c ~ 0.3 and similar initial eccentricity 
of the outer planet e c ~ 0.3 with simultaneously small eccentricity of the inner planet eb — 0. As we show 
in this paper, in such a case both systems would be located in a large island of the secular apsidal resonance 
with the semi-major axes antialigned in the exact resonance. There could be another planetary system, 
around HD 160691 (Jones et al. 2002), with a qualitatively similar dynamics as suggested by Gozdziewski 
et al. (2003). In a relatively small number of about 14 multi-planetary systems known to the date 4 , these 
three cases of supposedly identical configuration inspired us to ask whether such a dynamical resemblance 
can be casual or has a deeper origin in their formation scenario and/or orbital evolution. 

However recently, the Geneva team has published a new orbital solution substantially different than the 
first one (see Mayor et al. 2003). In the new solution the elements of the outer companion are dramatically 
altered. The period and the semi-major axis of the outer planet are changed by ~ 600 d and ~ 0.7 AU 
respectively. This obviously leads to a qualitatively different orbital configuration than the one analyzed 
before. The proximity of the HD 169830 and HD 12661 systems does not seem to be so clear anymore. Still, 
the orbital elements of the inner companions appear to be almost identical in the two planetary systems. 
Because the observational window of HD 169830 covers only about 1 orbital period of the outer planet, its 
orbital elements are still not well constrained. In our numerical experiments, we use both solutions although 
the old one essentially for a reference and an interesting comparison with the predictions given by the secular 
octupole- level theory (Lee & Peale 2003). 

According to Lee & Peale (2003), one has to employ a proper representation of the orbital fits to the 
RV observations, particularly in the case of the HPS. The orbital fits to the RV of HPS are best interpreted 
in Jacobi coordinates and these coordinates should be used in the dynamical investigations. Specifically, 
the authors demonstrate that the orbital time-evolution when expressed in the commonly used astrocentric 
elements can lead to spurious effects like a large scatter of the osculating semi-major axis and eccentricity 
of the outer planet about the secular value or a significant dependence of the integration results on the 
initial epoch. Unfortunately, the RV measurements of the HD 169830 system have not been published to 
the date. However, assuming that the best-fit parameters K,n 7 e,u),T p (where, for every planet involved, 
K is the semi-amplitude of RV variations, n is the mean motion, e is the eccentricity, u> is the periastron 
argument and T p is the time of periastron passage) describe a 2-Keplerian fit, they represent in fact the 
orbital parameters in the Jacobi coordinate system (Lee & Peale 2003) and using them, we can recover the 
semi-major axes a and the minimal masses of the companions. In some of the numerical experiments carried 
out in this paper, the inclination of the outer planet was changed by a tiny value, to 89.9°, to allow the system 
enter the third spatial dimension for the purpose of the numerical integrations. The periastron passage of 
the inner planet is selected as the initial epoch of the osculating elements inferred from the 2-Keplerian 
solutions. The two ICs used in the computations are given in Table 1. The first, preliminary one, IC is 
labeled by ICI and the current one by ICII. In some experiments we modified ICI, by setting e c ~ 10~ 6 , to 
avoid numerical singularities. The elements of the HD 12661 system are given in Table 3 and we refer to 
them as ICO throughout this work. 
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2. Numerical setup 

We use the so-called fast indicators, MEGNO and FMA as the main numerical tools to study the orbital 
stability of the HD 169830 system. The indicator called MEGNO (the Mean Exponential Growth Factor of 
Nearby Orbits), has been invented by Cincotta & Simo (2000) and we have applied it in the studies of the 
exosystems' dynamics described in a series of recent papers (for details and references on MEGNO see, for 
example, Gozdziewski et al. 2001; Gozdziewski & Maciejewski 2003). The MEGNO is a very efficient tool — 
typically, it makes it possible to distinguish between regular and chaotic dynamics during an integration of 
the system carried out over ~ 10 4 orbital periods of the outermost planet. This indicator provides a direct 
estimate of the stability in the strict sense of the maximal Lyapunov Characteristic Number (LCN). As in 
the previous works, the MEGNO integrations are driven by a Bulirsh-Stoer integrator, namely the ODEX 
code (Hairer & Wanner 1995). The relative and absolute accuracies of the integrator are set to 10~ 14 and 

5 • 10~ 16 , respectively. The zones of stability, corresponding to quasiperiodic motions of the system, are 
marked in the MEGNO maps by values close to 2. The positions of the nominal conditions are marked in 
the contour plots by the intersection of the two thin lines. In this work, we make yet another use of the 
integrations. During the MEGNO integrations we simultaneously analyze the osculating elements to study 
the short-term dynamics of the system. The osculating orbital elements feed the Frequency Map Analysis 
(FMA) algorithm (Laskar 1993). The FMA, an already classic fast indicator, allows us to detect orbital 
resonances and determine diffusion rates of the fundamental frequencies in the planetary system (Robutel 

6 Laskar 2001). The diffusion rate is a quantity that cannot be determined by the LCN, however it also 
measures the lack of regularity of an orbit. Along a quasi-periodic solution the fundamental frequencies 
are constant. For a chaotic orbit these frequencies change and by calculating their diffusion rates, one can 
directly detect the macroscopic changes of the orbital elements. Thus by combining both methods, it is 
possible to derive extensive information on the system's dynamics. 

Since our computations concern mostly a short-term dynamics of the planetary system, the FMA is 
very useful in identifying the positions of the mean motion resonances (MMRs). In order to identify these 
resonances, during every integration of MEGNO, the complex functions fk{U) = a*; expiAfc(ij) are computed 
for planets k — b,c at the discrete times ti, where tj+i — ij < 0.5Pb, over the time 2T of the order of 
10 4 P c (Pb, c are the orbital periods of the inner and outer companion, respectively). Here, ak denotes the 
semi-major axis and is the mean longitude of the planet. If the motion is quasiperiodic then in the time 
series {/fc(ij)}, the frequency Vk corresponding to the largest amplitude a° is one of the fundamental orbital 
frequencies of the system, called the proper mean motion, rik (Robutel & Laskar 2001). Because the secular 
frequencies are much smaller than the orbital frequencies, in the first approximation the precession of the 
orbits can be neglected and MMRs can be identified through the condition qv\> — pv c ~ 0, where p, q > 
are prime integers. The FMA gives the ratios v\,/v c — p/q, and then p and q can be easily resolved by the 
continued fraction algorithm. 

After finding v^l and v^fl where are ^ nc P r0 P er mean motions nb, c obtained after integrating 

the system over the time intervals [0, T], and [T, 2T], respectively, the diffusion rate is determined through 
(Tfe = 1 — v< k* l v k^ (Robutel & Laskar 2001). Following these authors, if the motion is quasiperiodic, Ok is 
equal to zero or, computed with the FMA algorithm, is a very small quantity, typically less than ~ 10~ 8 
while for an irregular chaotic motion, is larger by several orders of magnitude. The inferred diffusion of 
the semi-major axis , can be derived through the Kepler law n 2 a 3 — \i as An/n ~ — 3/2(Aa/a). In this 
way, although the time span of the integrations is relatively very short, the most relevant instabilities of the 
motion can be detected and identified. 
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In this paper, we use the FMA code kindly provided by David Nesvorny on his web page 5 . According to 
the author, the code incorporates the idea of the Frequency Modified Fourier Transform, FMFT (Sidlichovsky 
& Nesvorny 1997). We use a variant of the FMFT, with the so-called additional non-linear correction (see 
the cited paper for details). We have tested this code to make sure that it resolves the mean motions of 
non-interacting (Keplerian) orbits with the accuracy consistent with the internal precision of the method 
(set to 1(T 10 ). 

3. Orbital fits 

The two distinctively different orbital solutions, ICI and ICII, announced within a short period of time 
by the Geneva group strongly suggest that at least one of them corresponds to a local minimum of (xl) 1 ^ 2 - 
Unfortunately, it is not possible to verify their elements since the RV measurements are not available. In 
order to get at least some insight into the reliability of the newest set of orbital parameters, we have digitized 
the figure with RV measurements published on their WWW site 6 . The 93 data points we obtained differ 
from the real observations (for example, it is difficult to recover the exact moments of the observations), 
but they still properly describes the overall shape of the observed RV curve and its characteristic features. 
Note also that our set of "measurements" 7 is smaller than the real one (112) because in some parts of the 
RV curve it was very difficult to clearly resolve all points due to a dense sampling of the real measurements. 
We assume that the observational errors are given by a = <tq + G \ where er = 8.5 m/s and 0\ is an artificial 
normal random noise of mean and the standard deviation about 1, resulting in the "errors" in the range 
[6, 11] m/s and roughly approximating the scale of the real errors. 

We used the genetic algorithm scheme (GA) implemented by Charbonneau (1995) in his publicly avail- 
able 8 code PIKAIA to to look for a global, coplanar orbital solution. The data were modeled with 2-Keplerian 
model (as in Gozdziewski et al. 2003) and the self-consistent Newtonian model, i.e., driven by the full TV-body 
dynamics (Laughlin & Chambers 2001). For the first model the GA code was executed a few hundreds of 
times and repeatedly found a solution very close to the best one, given in Table 2. Even though the number 
of our "measurements" is smaller than in the real set, the similarity of the orbital solutions is very significant. 
Remarkably, the GA found very similar parameters for the outer planet (although the differences arc larger 
than for the inner planet). The best fit parameters have the (xt) 1 ^ 2 — 1-21 and the rms of ~ 10 m/s. This 
experiment was also repeated for two other data sets with 104 and 105 points (from two slightly different 
digitization) with the errors in the range [5, 12] m/s and [7, 10] m/s, respectively. The obtained best fit 
solutions were qualitatively the same as those given in Table 2 (let us note that the differences were most 
significant for e b ~ 0.29, e c ~ 0.29 and lo c ~ 270°). Our analysis indicates that the ICII by the Geneva group 
is really close to the global minimum of (xl) 1 ^ 2 an d should not change at least if based on the currently 
available set of the RVs. 

With the second model that accounts for the mutual interactions between planets, the best fit solution 
is not significantly better. It is characterized by the (xl) 1 ^ 2 — 1-21 and the rms of ~ 9.5 m/s. The 
corresponding orbital parameters are very similar to those from the 2-Keplerian fit (see Table 2; for a 

5 http: //www. boulder . swri . edu/ davidn/fmft/fmf t .html 
6 http : //obswww . unige . ch/ udry/planet/hdl69830_syst . html 

7 The set is available upon request from KG 
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reference, the synthetic RV-curves and the digitized "data" points are also shown in Fig. 4). This result 
suggests that the current RV set does not allow to detect significant TV-body coupling (for example, from 
low-order resonances). Finally, to roughly estimate the formal errors of the self-consistent Newtonian fit, 
we determined the confidence intervals of (xt) 1 ^ 2 (Press et al. 1992). The results are illustrated in Fig. 5. 
Clearly, the formal errors of a c and e c (semi-major axis and eccentricity of the outer planet) are substantial 
but the minimum of (xt) 1 ^ 2 is weu localized. We note that to estimate the conservative errors in the best-fit 
parameters, the effect of stellar "jitter" should be taken into account (as in Gozdziewski & Maciejewski 
2003). 

4. Orbital stability of the HD 169830 system 

Both orbital configurations of HD 169830 , ICI and ICII, appear to be stable. The evolution of the 
Jacobi osculating orbital elements computed over the time span of ~ 3.6 • 10 4 revolutions of the outer planet 
(about 150,000 yr) for the ICI fit is illustrated in Fig. 1. During this time, the eccentricities (see Fig. lc) 
vary with a large full-amplitude of about 0.35. Fig. Id reveals a presence of the secular apsidal resonance 
(SAR) in which the apsides are on the average antialigncd. The semi-amplitude of the librations, 6 (where 
6 = tub — vj c and w\, tC are the respective longitudes of periastron), is about 90°. Apparently, the orbital 
elements vary in a regular way and this is strictly confirmed by the MEGNO signatures. The temporal 
changes of MEGNO, Y(i), are shown in Fig. le and its mean value, (Y)(t) (Fig. If), perfectly converges to 2. 
For a comparison, we computed these quantities for the HD 12661 system using the ICO from Table 3. The 
results are shown in Figure 3. In both these cases, the variations of the orbital elements and the character 
of MEGNO convergence are qualitatively identical. However, the ICII fit for the HD 169830 system results 
in a qualitatively different orbital evolution (Fig. 2). The SAR docs not seems to be present anymore. 
Nevertheless, the MEGNO signature (computed over about 550,000 yr) indicates a stable, quasiperiodic 
configuration. 

The MEGNO signatures have been verified by direct 1 Gyr integrations using the RMVS3 integrator 
from the SWIFT package (Levison & Duncan 1994). We repeated the integrations with two different time 
steps, equal to 8 and 10 days. As one would expect, no instability occurs during this time and the orbital 
elements vary within the bounds determined by the short-term integrations. We do not show these results as 
they are basically an extension of the plots shown in Figures 1, 2 and 3. Obviously, due to the uncertainties 
of the orbital fits, such an examination of the isolated IC's is not representative for the system's dynamics. 
In order to find out whether the dynamics is robust to small adjustments of the initial condition, we have 
computed one-dimensional scans of (Y) by changing the semi-major axis of the outer planet and keeping the 
other orbital parameters fixed at their initial values given in Table 1 for HD 169830 and 3 for HD 12661 . 
The results of this experiment are shown in Fig. 6. The MEGNO scan (the upper graphs of the respective 
panels in Fig. 6), computed with the resolution of about 3 • 10~ 4 AU, reveals a number of spikes. Most of 
these spikes represent the MMRs between the planets. This identification is based on the FMA as described 
in the previous Section. We marked MMRs of the order p + q not grater than about 20. These scans reveal 
that the dynamical environment of the nominal ICs is qualitatively the same for the ICI of the HD 169830 
and the ICO of HD 12661 systems. The outer planet lies in a stable zone between strong MMRs 6:1 and 7:1. 
The graph for the ICII of the HD 169830 system reveals a different dynamical environment. Clearly, the 
change in e c (compared to e c from ICI) "shifted" the system to the zone between the 9:1 and 10:1 MMRs. 

The FMA working with the MEGNO code enabled us to find out how sensitive both methods are to 
the instabilities of the motion. Bottom parts of the panels of Fig. 6 show the diffusion rate computed for 
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the outer planet and estimated over its ~ 18, 000 orbital periods. In the regions of quasiperiodic motions, as 
classified by MEGNO, the diffusion rate, cr c , is smaller than about 10~ 8 . Such small values ensure us that 
the motion is close to a quasiperiodic evolution. In the MMRs zones a c grows up to 10~ 2 — 10 _1 . Hence in 
these areas, the osculating orbital elements exhibit macroscopic changes. Both algorithms are in an excellent 
accord, as they provide the same positions of the MMRS and very similar estimates of their widths. In other 
experiments, we noticed that for much shorter integration times, ~ 0.45 • 10 4 P C , both algorithms can still 
detect all the relevant MMRs although in this case they do not allow to point out clearly enough some of 
the weak resonances. In fact, for the four times longer integration (~ 1.8 • 10 4 P C ), the FMA seems to be 
even more sensitive to the presence of weak resonances than the MEGNO algorithm is. We note here that 
by a compromise forced by the numerical efficiency of the code, in further runs of MEGNO, we have set the 
integration time to about 0.9 • 10 4 P C . 



5. Global dynamics of the HD 169830 system 

To extend the above one-dimensional analysis, we can investigate the stability of the HD 169830 system 
in a few representative planes of its orbital parameters. 

The program computing MEGNO simultaneously evaluated the maximal values of the eccentricities, 
ej^ x , attained during the integration time. We also stored the maximal value of semi-amplitude of the 
librations, # max 7 after every step (set to P c ) of renormalization of the variational equations. It helped us to 
detect the apsidal resonance and to estimate the semi-amplitude of the librations. The maximal value of 
the critical argument 9 was taken relative to the center of the libration 0° or 180°. To avoid the effects of a 
possible transition into the SAR, the determination of # max was started after the first half of the integration 
period (about 0.45 • 10 4 P C ). Finally, if 8 max < 90°, then we treated this value as a semi-amplitude of the 
apsidal librations. The period of the integrations is relatively short but as we show in Gozdziewski (2003), 
such information on the short-term dynamics can still give us much insight into the global behaviour of the 
system. The results are illustrated in Figs. 7-9, where the left panels are for MEGNO, (Y), the middle 
panels are for # max and the right panels are for e™ ax . In these scans, the initial parameters that are varied 
are the coordinates of the maps and the other initial orbital elements are fixed at their nominal values given 
in Table 1. 

Figure 7 is for the (a c , e c )-plane. Since the resolution of these maps is 200 x 50 data points and the 
integration time is shorter than that of the one-dimensional a c scan (shown in Fig. 6), there is a lack of 
some fine resonance structures visible in Fig. 6. Yet the location of the dominant low-order MMRs: 6:1, 
13:2 and 7:1 for ICI and 8:1, 9:1 and 10:1 for ICII, as well as their widths are clearly marked. The system 
would be chaotic for e c roughly grater than 0.5-0.6. In the central, stable parts of the MMRs, the maximal 
values of eb are small indicating their stabilizing influence on the motion. For ICI of the system around 
HD 169830 (as well as for ICO of the HD 12661 system, not shown here), the plot for # max reveals an extended 
zone of the SAR about the libration center of 180°. For very small initial e c ~ 0, the semi-amplitude of 
librations reaches the limiting 90°. If the eccentricity becomes larger than ~ 0.005 then the semi-amplitude 
decreases rapidly to 60°-70°. The same type of (a c , e c )-scans of MEGNO for the HD 12661 system is shown 
in Gozdziewski (2003) (his Fig. 6) and Gozdziewski & Maciejewski (2003) (their Fig. 3). In both cases, the 
MEGNO structures around the nominal ICs are very similar. Interestingly, the (a c , e c )-maps of MEGNO for 
some of the fits to the RV data of the HD 160691 planetary system in Gozdziewski et al. (2003) (e.g., the 
fits GM5 and GM6, Fig. 5 in that paper), describe a dynamical setup very similar to these of the HD 169830 
(with ICI) and HD 12661 systems. This is no longer true for the ICII of the HD 169830 system. In this 
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case, the zone of the SAR is confined to a very small e c (but over wide range of a c ) and to the centers of the 
MMRs (see Fig. 7, bottom panels). 

Figure 8 shows the behaviour of MEGNO, the semi- amplitude of apsidal libration and the maximal 
eccentricity of the inner planet for HD 169830 system on the (eb, e c )-planc. The nominal configuration of 
the system (denoted as previously by the intersection of the two thin lines) is located in an extended zone 
of quasiperiodic motions as seen in the MEGNO map (Fig. 8, panels in the left column). The e™ ax -maps 
(Fig. 8, right column) enable us to determine the border of the global instability established by these values 
of the initial e c which lead to large maximal eccentricities. For ICI, this border is substantially shifted 
towards larger initial eccentricities compared to the border of chaotic zone visible in the MEGNO map but 
overall it has a similar shape. The same effects are seen for ICII. In the SAR-maps (Fig. 8, middle column), 
for both IC's one can find a very sharp border of this resonance in the region of small eb- For ICI, it 
ends at the border of unstable motion present in the relevant e™ ax panel. For ICII, the semi-amplitude of 
the librations is about 60°-70° in this area but as for ICI (Fig. 7, middle panel), a very narrow strip of 
semi- amplitudes of ~ 90° for the initial e c ~ is present. These large- amplitude librations for small e c are 
explained in the next section. Let us note that also in this test, the obtained picture of the stability zones 
for ICI of the HD 169830 system is qualitatively the same as for the HD 12661 system — see Fig. 5 and 
the discussion from Gozdziewski (2003). However, the SAR- map is quite different for ICII since the SAR is 
confined to a narrow zone of regular motion about small initial e c and also to unstable zone which is seen 
around (e b ,e c ) = (0.47,0.47). 

In the next set of experiments, we determine the zones of stability while the system inclination and 
hence planetary masses are varied. Assuming that the (unknown) inclinations of the orbits are the same 
i = ?b = ic, the relative inclination i le \ of the orbits depends on the difference of the nodal longitudes, 
AO = O c — Ob, according to the formula: 

cos(i TO i) = cos(ib) cos(z c ) + sin(ib) sin(z c ) cos AO. 

We assume that both orbits are prograde. By changing i, the masses are altered because the minimum 
masses, mb, c sini, corresponding to edge-on orbits, must be preserved. 

The results are shown in Fig. 9. For both ICI and ICII, even for large relative inclinations, the zones 
of stability extend up to very low system inclinations, i ~ 15°. This is supported by the ej^^-scan which is 
shown in the right panel of this figure. In some unstable zones, seen in the MEGNO map for the ICI (Fig. 9, 
top-left panel), we can identify the 13:2 MMR. It quickly destabilizes the system in the areas about i ~ 45° — 
in these regions both eccentricities increases rapidly up to 0.8-0.9. In the strip about low inclinations, most 
likely thanks to the SAR, this effect is absent. 

For ICI, due to the initial e c ~ and a geometrical singularity (for very small eccentricity the argument 
of periastron becomes undefined) , our numerically computed # max -map (top- middle panel of Fig. 9) does not 
allow to clearly detect the zone of the SAR. However, in our preliminary study of the ICI we set e c ~ 0.007 
and then we found a wide zone of the SAR in the central part of the scanned area. The overall picture is 
qualitatively the same as in the case of the HD 12661 system (Gozdziewski & Maciejewski 2003) where for 
low-z re i configurations, the SAR persists almost in the entire range of the system inclination i. 

For ICII, the unstable zones in the MEGNO map are basically similar to those of the ICI of the 
HD 169830 and HD 12661 systems. Interestingly, all the maps for maximal eb reveal characteristic quarter- 
circle areas in which the eccentricity grows up to 1 . In the strips near the central parts and corresponding to 
eb — 1 the configuration becomes chaotic (compare the MEGNO- and ej^-maps) . The analysis of the orbits 
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in these areas shows that together with the excitation of the inner planet's eccentricity, the inclination of this 
planet oscillates with a very large amplitude. It is illustrated in Figure 10 for the ICII and the initial relative 
inclination set to 80°. Let us observe that in this case e b approaches 0.9 while the inclination goes up to 
140°. Nevertheless, the system still remains rigorously stable. This phenomenon seems to be an analogue of 
the Kozai resonance (Kozai 1962) that is known for highly inclined asteroids and generally has been observed 
in triple systems in which the mass of the inner body is much smaller than the mass of the central and the 
outer body (Holman et al. 1997). In the multi-planetary exosystems, the masses of the perturber and the 
inner body are comparable which makes the problem more complex than in the asteroidal approximation. 
Let us note that a possibility of this phenomenon has been analyzed in v Andromedae case (Chiang et al. 
2001). The inclined configurations which we study here are rather special because it has been assumed that 
the initial inclinations are the same for both planets. However, a similar behaviour can be expected when 
the absolute inclinations are varied independently. It has been observed also in other hierarchical exosystems 
(see for example Gozdziewski 2003; Gozdziewski 2003). Since this phenomenon provides strong dynamical 
limits on the relative inclinations, it certainly deserves a detailed study which we plan to carry out in a 
future paper. 



6. Secular apsidal resonance 

The SAR in the HD 12661 exosystem was analyzed numerically in Gozdziewski (2003). We found an 
extended zone of the SAR in wide ranges of the semi-major axes, eccentricities and masses of the companions. 
Similar behavior of the SAR can be observed for the ICI of the HD 169830 system. It should be noted that 
those results were based on the IC of the HD 12661 system derived from the 2-Keplerian model by Fischer 
et al. (2003). Our attempts to analyze the same RV data set resulted in a rather different orbital solution 
(Gozdziewski & Maciejewski 2003). Yet, the overall SAR features did not change. These results have been 
confirmed and greatly extended analytically by Lee & Peale (2003) in the framework of their octopole-level 
secular planetary theory. 

Here, we literally quote that basics of their theory, which are relevant for our further analysis. This 
theory describes the secular dynamics in a coplanar hierarchical 2-planet system. It is assumed that MMRs 
are absent. The equations of motion originate from the Hamiltonian which is expanded up to the third order 
in parameter a = ab/a c and averaged over both mean longitudes lj, j = b,c. The relevant parameters of 
the secular theory are defined through the momenta conjugated to the mean longitudes and arguments of 
periapse gj — u)j : 

L b = ■ \/G{m* +m h )a b , 

ra* + m b 

(m* + m b )m c 



L c = ; ■ y G(m± + m b + m c )a c , (1) 

+ mb + m c 



G J =L j jl-ej, 

where to* is the mass of the parent star. The momenta are, respectively, the magnitude of the maximum 
possible angular momentum (for circular orbits) and the magnitude of the angular momentum. The sum of 
the remaining two momenta, Hj = Gj cos ij , is the z-component of the angular momentum of the system 
conjugate to the longitudes of the ascending node hj = Clj. Let us note, that this formulation is valid in 
the general case of an inclined system. In coplanar configurations u>j are undefined but then the longitudes 
of periastron Wj take over their role. In the averaged system, lj are absent in the Hamiltonian and Lj are 
constants of motion. Because the octopole-level Hamiltonian depends only on the combination w b — zu c , it 
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turns out that the total angular momentum Gb + G c is also the integral of motion in the averaged system. 
This reduces the coplanar dynamics to 1 degree of freedom, with eb or e c , and 9 as the relevant phase-space 
variables. The authors found that for the same constants: 

_ 5 (ro* - m b ) a 
4 (m* + m b ) 

A = Lb/L c , the same initial eb,e c and 9 = — m c the averaged equations of motion describe the same 
trajectories in the phase-space diagram of eb jC (#)- For planetary masses much smaller than the central body, 
these trajectories and the amplitudes of eccentricity oscillations should be independent of the inclination 
of the system. Further, by introducing the parameter 7 = (Gb + G c )/(Lb + L c ) which describes the non- 
dimensional total angular momentum of the system, one defines the critical value of A, 

crit ~ 5 - 3 7 2 ' 

According to the equations of motion in the averaged system, if A ~ A cr it then libration of 9 is almost certain 
with possibly large amplitude variations of both eccentricities. In general, the librations can take place about 
two libration centers, 0° or 180°. Using the fits published in Fischer et al. (2003), the authors found that 
indeed, for HD 12661 system, A ~ A cr it and they identified an extended libration island of the SAR with the 
apsides antialigncd in the exact resonance. 

Following Lee & Peale (2003), we examined the ICs for the HD 169830 and HD 12661 systems by 
calculating these phase-space diagrams numerically. To obtain such diagrams, we fixed the constant level 
of the total angular momentum corresponding to the nominal IC and the phase curves were computed for 
varied e c (eb was calculated from the total momentum integral, the formula 1). The time of integrations 
was the same as for the derivation of the 2-dimensional stability maps. The results are shown in Figure 11. 
For both ICI of the HD 169830 and ICO of the HD 12661 systems, the phase space is occupied by two large 
libration islands about 9 = 0° and 9 = 180° which is consistent with the secular theory. For the ICI of the 
HD 169830 system, A = 0.644 and A cr i t = 0.9, so these parameters substantially differ. However, Fig. 11 
shows that the occurrence of a libration mode still appears more likely than the 9 rotations. For the ICO of 
the HD 12661 , A = 0.75 is closer to its critical value A cr i t = 0.87, and indeed, the the resonance island is 
more extended. 

The configuration corresponding to ICI of the HD 169830 system is localized almost on the border of 
the island of the antialigncd SAR. This border is a separatrix and if the initial e c ~ 0, the semi- amplitude 
of 9 librations approaches 90° as the system enters the vicinity of the separatrix. In the resonance island, 
the librations are about the libration center 9 — 180°. The diagram also helps us to understand the large 
amplitudes of the eccentricities. Qualitatively, one obtains the same picture for the HD 12661 system. For 
the ICII of the HD 169830 system, we have A = 0.97 and A crit = 0.77. In this case the nominal HD 169830 
system lies in an extended area of ^-rotations seen in Fig. 11. 

The octupole theory allows us to quickly examine the extent of the SAR in the space of the orbital pa- 
rameters. Because the elements of the inner planet are well determined, the overall picture of the HD 169830 
dynamics depends mostly on the not too well constrained parameters of the outer planet. Assuming that 
strong mean motion resonances are absent, the secular dynamics can be investigated in detail through solving 
Equations (1) — (4) of Lee & Peale (2003). These equations describe the secular time-evolution of eb, c and 
H7b iC - Unfortunately, there seems to be no explicit, analytical solution to these equations. Also a construc- 
tion of a simple SAR criterion, like the one developed by Laughlin et al. (2002) or Ji et al. (2003) in the 
framework of the Laplace-Lagrange secular theory, does not appear to be possible, cither. However, these 



- 10 - 



equations are very simple for a numerical treatment and their direct integration is rapid, hundreds times 
faster than the TV-body integrations. The SAR can be detected in the parameter space by looking for the 
maximal deviations of 8 from the given libration center 0° or 180°. It is the same method of detecting the 
SAR as the one applied in the full A-body integrations but here it is well justified because according to the 
generic secular dynamics, the SAR appears about two libration centers. Assuming that the integration time 
covers at least one secular period (so as it is not too short), we can easily detect the librations or rotations 
of the critical argument 9. 

Using this approach, we computed a number of the SAR-maps in the (a c ,e c ) and (eb,e c ) domains 
for the initial conditions ICI, ICII and ICO. In order to directly compare the results, the ranges of the 
parameters were the same as for the numerical maps. The results are shown in Fig. 12. Generally, they 
agree in the regions corresponding to small eccentricities, for which both the shape of these regions and the 
semi-amplitude of librations are well reproduced by the secular theory. However, we can also find significant 
discrepancies between the analytical and numerical estimations. For example for ICII, in the (a c , e c )-planc 
one can find islands of the SAR which are absent in the map computed by solving the secular equations. 
The differences are also apparent in the (eb, e c ) plane for all examined IC's. In these cases, the results agree 
only in the regions of relatively small e c . 

In order to understand these differences, we compared the (e c ,6>) diagrams obtained by the direct 
numerical integrations and by solving the secular equations. The tests were performed for: (a) ICI in which 
we changed eccentricities to eb = 0.3, eb = 0.3, (b) ICII with a c = 3.51, e c = 0.47 (see bottom-middle panel 
in Fig. 7 corresponding to the center of the 9:1 MMR) (c) ICII with (eb = 0.47, eb = 0.47) (corresponding 
to the island in the bottom middle panel in Fig. 8). The results are presented in Fig. 13. In all these cases, 
the full dynamics seems to be only roughly described by the secular approximation. The disagreement is 
obviously expected in the case (b) since the assumptions of the secular theory are violated (an MMR is 
present). In other two examples, the differences are most likely due to significant mutual interactions caused 
by large eccentricities resulting in deformations of the phase-space curves (case a) or a chaotic evolution 
(case c). Clearly the secular theory should be applied with caution. 

7. Conclusions 

In this paper we carry out a dynamical analysis of the recently discovered 2-planetary exosystem around 
HD 12661 . According to the preliminary initial orbital parameters announced by the Geneva Extrasolar 
Planet Search team (orbital solution ICI), this system resembled the other known planetary hierarchical 
exosystem HD 12661 and the HD 160691. The recently updated orbital solution (ICII) has been changed 
in a significant way and the close similarity of the HD 169830 and HD 12661 systems is not so apparent 
anymore. 

The exosystems around HD 122661, HD 160691 and HD 169830 belong to a class of the hierarchical 
planetary systems. They are characterized by a relatively small ratio of the semi-major axes, ~ 0.1 — 0.3. 
Their dynamics should be analyzed in the Jacobi orbital elements (Lee & Peale 2003; Gozdzicwski et al. 
2003). The commonly used astrocentric elements may lead to small shifts of the positions of the orbital 
resonances and to differences in the evolution of the osculating orbital elements. Following these papers 
we used the two orbital fits published by the discovery team but we recomputed the semi-major axes and 
masses of the companions. The current ICII was verified by our independent analysis of the approximate 
RV measurements which we obtained by digitizing a figure with the synthetic best-ft RV curve and the real 
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observational points. Such an "approach" to obtain the "observations" was forced by the lack of access to 
the real data. Since the (xl) 1 ^ 2 function will typically have many local minima, we cannot be sure that a 
solution is the proper one without a global (xl) 1 ^ 2 analysis. As demonstrated on the solution ICI and ICII 
of the HD 169830 system, the difference between the global and a local minimum may lead to enormous 
differences in the overall dynamical picture of the system. Obviously, the digitization cannot provide very 
accurate moments of the observations, nevertheless such "data" describe the the overall shape of the real 
RV measurements quite precisely and make it possible to obtain some insight into the quality of the best 
fit solution by the discovery team. Using the genetic algorithm to find the global minimum of {x 2 ) 1 ^ 2 an d 
describing the measurements with the 2-planet Keplerian and Newtonian models, we found almost the same 
orbital elements of the inner planet as these reported by the discovery team and very similar elements for the 
outer companion. The application of the 2-Keplerian and TV-body models produced the same value of (xt) 1 ^ 2 
and similar orbital solutions, ft demonstrates the absence of strong interactions between the planets, ft also 
favors the solution with a large separation between the companions rather than the previously announced 
configuration. Apparently, the first orbital solution was just a local minimum of the (xt) 1 ^ 2 function. 

The two orbital solutions have different dynamical features. The ICI solution is close to an unstable 
13:2 MMR and lies between 6:1 and 7:1 MMRs. The exosystem corresponding to ICII is located between 
9:1 and 10:1 MMRs. Both configurations appear to be stable on the Gyr time scale as demonstrated by 
our long-term integrations, the fast indicator analysis and the secular octupole-level theory. In the phase 
space, both ICs are lie in the wide regions of stable motions in spite of large variations of the eccentricities 
(up to 0.5). Obviously, these conclusions may changed in the future as the time-span of the current RV 
measurements is shorter than the period of the outer planet. Moreover, the moderately constrained orbital 
elements of the outer planet do not allow us to exclude a proximity of the system to an unstable low-order 
resonance. The main dynamical difference between the two ICs is the SAR with the apsidal lines antialigncd 
in the exact resonance. It is present for ICI in the wide ranges of the orbital parameters while for the ICII 
it seems to be ruled out. In the first case, the HD 169830 would be very similar to the hierarchical resonant 
systems. In this "class" we can find the discussed HD 12661 and possibly the HD 160691 exosystem. The 
HD 37124 (Butler et al. 2003) system could also be a member of this group. However, in the HD 37124 system 
the apsides are most likely aligned (Ji et al. 2003). If the two outer planets are taken into consideration, the 
v Andromedae (Butler et al. 1999) system could also be assigned to this group (Chiang & Murray 2002). 
However, the new orbital solution ICII favors the proximity of the HD 169830 system to the HPS with the 
outer, very massive (possibly sub-stellar) companion like the HD 38529 (Fischer et al. 2003), HD 74156 9 
and HD 168443 (Marcy et al. 2001; Udry et al. 2002) systems. In these systems, the stability is maintained 
thanks to a relatively large separation of the companions. Nevertheless, due to very large masses involved 
and large amplitudes of the eccentricity variations, the mechanism providing the stability is puzzling. 

The results of the numerical analysis of the orbital evolution and stability can be compared with the 
predictions by the secular octupole-level theory of Lee & Peale (2003). The results of both approaches are in 
accord for moderately low eccentricities of the outer planet. The discrepancies appear if both eccentricities 
are large in the regions of the phase space close to the MMRs and in the chaotic zones. It seems that 
although the secular theory makes it possible to explain the main features of the HD 169830 -like systems, 
the direct integrations are still necessary to understand the dynamics in detail. 

For a better grasp on the dynamical picture of the HD 169830 system, we need to wait for some time, 
ideally until the observations cover a few orbital periods of the outer companion. Since the orbital elements 



9 http : //obswww.unige . ch/ udry/planet/hd74156 .html 
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of the inner companions are similar in the HD 169830 and HD 12661 systems and are already determined 
with high accuracy, there is still a chance that our preliminary hypothesis that the two planetary systems 
are dynamically related is plausible. 
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Fig. 1. — The orbital evolution (expressed in the osculating orbital elements of Jacobi) of the HD 169830 
system for the nominal ICI given in Table 1. Panels (a,b) show changes of the semi- major axes. Panel (c) 
is for the eccentricities. Panel (d) shows the critical argument of the apsidal resonance 9 = — vj c . Panel 
(e) is for MEGNO as a function of time, Y(t), and panel (f) is for its mean value (Y). 

Fig. 2. — The orbital evolution of the HD 169830 planetary system for the ICII given in Table 1. The 
description of the panels is the same as in the previous figure. 

Fig. 3. — The orbital evolution of the HD12661 planetary system for the IC given in Table 3. The description 
of the panels is the same as in the previous figure. 

Fig. 4. — The digitized RV data of the HD 169830 system and the synthetic RV curves corresponding to the 
best fits found from this paper. The thicker line is for the TV-body model. The orbital parameters of the 
solutions are given in Table 2. 

Fig. 5. — Confidence intervals corresponding to la, 2a and 3a levels of (xt) 1 ^ 2 obtained for the best TV-body 
fit to the digitized RV data. The fit is given in Table 2. 

Fig. 6. — Scans of MEGNO (upper graphs in every panel) and the diffusion rate of the proper mean motion 
of the outer planet (lower graphs labeled by loga c ) for the ICI of the HD 169830 system (the top panel) 
ICII of the HD 169830 system (the middle panel) and for the HD 12661 system (the bottom panel) over the 
semi-major axis of the outer planet (in AU). The resolution of the plot is ~ 3 • 10~ 4 AU. Labels mark the 
positions of the strongest mean motion resonances (up to the order about 20) as they have been identified 
by the FMA. The nominal values of a c are marked with big dots. 

Fig. 7.— The stability of the HD 169830 system in the (a c , e c )-planc. The left panels are for MEGNO. 
The middle panels give an estimate of the semi-amplitude of the apsidal librations, 6» max ; about the line of 
apsides antialignment (note that the white areas are for 6* max > 90°). The right panels show the maximal 
eccentricity of the inner planet. The data grid has the resolution of 0.006AU x 0.005. The top row is for the 
ICI, the bottom row is for the ICII. 

Fig. 8. — The stability of the HD 169830 system in the plane of the eccentricities. The left panels are for 
MEGNO. The middle panels give an estimate of the semi- amplitude of the apsidal librations, # max , about 
180°. Note that the white areas are for 6* max > 90°. The right panels show the maximal eccentricity of the 
inner planet. The data grid has the resolution of 100 x 100 points. The top row is for the ICI, the bottom 
row is for the ICII. 

Fig. 9. — The stability map of the HD 169830 system when the relative inclination of the orbits and 
the absolute inclination i are varied. The left panels are for MEGNO. The middle panels show the semi- 
amplitude of librations of 9 about 180°. The white areas are for # max > 90°. The right panels show the 
maximal eccentricity of the inner planet. The data grid has the resolution of 3° x 1°. The top row is for the 
ICI, the bottom row is for the ICII. 

Fig. 10. — The evolution of the inclination and eccentricities in the planetary configuration corresponding 
to ICII when the initial relative inclination of orbits has been changed to 80° . The left panel is for the inner 
planet and the right panel is for the outer companion. The bottom plots are for the eccentricities magnified 
by 100. 
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Fig. 11. — The phase-space diagrams at the same level of the total angular momentum. The left column is 
for the ICI, the middle column is is for the HD 12661 system, the right column is for the ICII. The nominal 
ICs are marked with a big dot. The phase-space trajectories have been computed for the initial 9 = 0° 
and for 9 = 180°; e c varies while eb is determined from the integral of the total angular momentum. Other 
orbital parameters are fixed at their nominal values. 

Fig. 12. — The semi-amplitude of the SAR about the libration center 180° in the planes of semi-major axes 
and eccentricities as predicted by the secular octupole theory. The left column is for the ICI, the middle 
column is for the ICII of the HD 169830 system and the right column is for the ICO of the HD 12661 system. 
The white areas are for fl max > 90°. 

Fig. 13. — The comparison of the phase-space diagrams at the same level of the total angular momentum 
obtained from the numerical integration of the full equations of motion (top row) and derived from the 
secular-octupole theory (bottom row). The initial conditions correspond to the modified ICI and ICII (see 
text for explanation). The modified initial conditions are marked with big dots. They correspond to the 
coordinates marked with (a), (b) and (c) in Fig. 7, 8, respectively. 
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Table 1: Jacobi orbital parameters of the HD 169830 system from 
http://obswww.unige.ch/udry/planet/hdl69830_syst.html (updated on August, 28, 2003; the 
ICII solution). The first version of the fit (dated 28 June, 2003; the ICI solution) is given in the parentheses. 



The mass of the central star is 


equal to 1.4 M©. 




Jacobi orbital parameter 


Planet b 


Planet c 


ni2 sini [Mj] 


2.88 (3.03) 


4.05 (2.51) 


a [AU] 


0.811 (0.816) 


2.856 (3.598) 


P [d] 


225.62 ±0.22 (227.43) 


2102 ±264 (1487) 


e 


0.31 ±0.01 (0.327) 


0.33 ±0.02 (0.0) 


w [deg] 


148 ±2 (156.1) 


252 ±8 (44.0) 


T p [JD-2,400,000] 


51923 ±1 (51472.43) 


52516 ±25 (50101) 


M(T b ) [deg] 


0.00 


258.4 (332.2) 



Table 2: Jacobi orbital parameters of the HD 169830 system based on the digitized figure from 
http://obswww.unige.ch/ udry/planet/hdl69830_syst .html (the version updated on August, 28, 2003). 
Numbers without parentheses are for Jacobi elements of the 2-Keplerian solution. Values in parentheses are 
for the osculating, Jacobi elements from the iV-body, self-consistent fit given for the date of the first digitized 
observation (JD=2,451,294.97). The mass of the central star is equal to 1.4 M . 



Orbital parameter Planet b Planet c 

m 2 sini [Mj] 2.86 (2.85) 3.88 (3.84) 

a [AU] 0.812 (0.812) 3.41 (3.38) 

P [d] 225.62 1939.4 

e 0.307 (0.307) 0.334 (0.336) 

a; [deg] 148.6 (148.8) 260.1 (260.5) 

T p [JD-2,400,000] 51922.5 52533 

M [deg] 0.0 (78.9) 204.7 (128.1) 

(Xl) 1/2 1-21 (1.21) 

RMS [m/s] 10 (9.5) 

Vp [m/s] 3.35 (6.3) 



Table 3: Orbital parameters of the HD 12661 system from Gozdziewski & Maciejewski (2003) for the epoch 
of the first observation (JD=2,450,831.608). The mass of the central star is equal to 1.08 M©. The orbital 
periods come from the FMA as reciprocals of the proper mean motion frequencies. 



Orbital parameter 


Planet b 


Planet c 


ni2 smi [Mj] 


2.33 


1.69 


a [AU] 


0.82 


2.78 


P [d] 


263 


1632 


e 


0.349 


0.076 


w [deg] 


115.2 


294.4 


M [deg] 


129.37 


352.9 



